##########################################
## "Voters get what they want"          ##
##########################################
## Heinrich, Kobayashi, Long            ##
##########################################

##          
## Descriptives of Nielsen data
## Script 09
## June 29, 2017


## Load data
############
load("output/data_SH.Rdata")
load("output/data_EH.Rdata")


## Specifications for subsets
#############################
A <- c("Leader news", "HR news")
B <- c("Economic hierarchy", "Security hierarchy")
C <- data.frame(Name=c("All data", "Above 40th p'tile", "Below 40th p'tile"), 
                Lower=c(0.0, 0.6, 0.0),
                Upper=c(1.0, 1.0, 0.4))


## Calculate correlations
#########################
corrs <- c()
for(i in 1:length(A))
{
  for(j in 1:length(B))
  {
    for(k in 1:nrow(C))
    {
      if(B[j] == "Economic hierarchy") this_data <- data_EH
      if(B[j] == "Security hierarchy") this_data <- data_SH
      
      if(A[i] == "Leader news")
      {
        this_data <- subset(this_data, llnNYT_leader >= quantile(this_data$llnNYT_leader, C$Lower[k]) &
                              llnNYT_leader <= quantile(this_data$llnNYT_leader, C$Upper[k]))
      }
      if(A[i] != "Leader news")
      {
        this_data <- subset(this_data, llnNYT_HR >= quantile(this_data$llnNYT_HR, C$Lower[k]) &
                              llnNYT_HR <= quantile(this_data$llnNYT_HR, C$Upper[k]))
      }
      
      c_hier_HR <- cor(this_data[, c(ifelse(B[j] == "Economic hierarchy", "us_EH1995", "us_SH1995"), "lphysint")])[1, 2]
      c_news_HR <- cor(this_data[, c(ifelse(A[i] == "Leader news", "llnNYT_leader", "llnNYT_HR"), "lphysint")])[1, 2]
      c_news_hier <- cor(this_data[, c(ifelse(A[i] =="Leader news", "llnNYT_leader", "llnNYT_HR"), 
                                       ifelse(B[j] == "Economic hierarchy", "us_EH1995", "us_SH1995"))])[1, 2]
      c_news_news <- cor(this_data[, c("llnNYT_leader", "llnNYT_HR")])[1, 2]
      
      corrs <- rbind(corrs, 
                     data.frame(NewsType=A[i],
                                Hier=B[j],
                                NewsSet=C$Name[k],
                                Correlation=c("Cor(Hierarchy, Violations)",
                                              "Cor(News, Violations)",
                                              "Cor(Hierarchy, News)",
                                              "Cor(News/NYT, News/Human rights)"),
                                value=c(c_hier_HR, c_news_HR, c_news_hier, c_news_news)))
      
    }
  }
}


## Make graph
#############
corrs$News <- paste0(corrs$NewsType, "\n", corrs$NewsSet)
corrs$Hier <- str_replace_all(string=corrs$Hier, pattern=" hierarchy", replacement="")
g <- ggplot(data=corrs, aes(x=Correlation, y=value, shape=Hier))
g <- g + geom_hline(yintercept=0, size=.6) + xlab("")
g <- g + geom_point(size=2.7, alpha=.7, colour="black", fill="black", position=position_dodge(width=.24))
g <- g + facet_grid(. ~ News)
g <- g + scale_shape_manual(values=c(21, 23), name="Hierarchy")
g <- g + theme_bw() + ylab("Correlation") + coord_flip()
g <- g + ggtitle("Correlations of key variables by samples")
g <- g + theme(axis.text = element_text(size=rel(.59)),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),               
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               strip.text = element_text(size=rel(.74), hjust=0),
               strip.background = element_blank(),
               plot.title = element_text(size=rel(0.89), face="bold", hjust=0))
g
ggsave(file="output/figures/A1-Descriptives.pdf", width=10, height=6)

    





    
    
